{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation and Differentiation\n",
    "\n",
    "All algorithms use the baryentric weights as described in the book \"Implementing Spectral Methods for PDEs\" by David Kopriva.\n",
    "\n",
    "- Interpolation can be performed either with `interpolate(dest, u, basis)` or via an \n",
    "  `interpolation_matrix(dest, basis)`.\n",
    "- In a nodal basis, the derivative matrix is available as `basis.D`. Alternatively, `derivative_at(x, u, basis)`\n",
    "  can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "# define nodal bases\n",
    "p = 5 # polynomial degree\n",
    "basis1 = LobattoLegendre(p)\n",
    "basis2 = GaussLegendre(p)\n",
    "\n",
    "# the function that will be interpolated\n",
    "ufunc(x) = sinpi(x); uprim(x) = π*cospi(x)\n",
    "#ufunc(x) = 1 / (1 + 25x^2); uprim(x) = -ufunc(x)^2*50x\n",
    "\n",
    "for basis in (basis1, basis2)\n",
    "    u = ufunc.(basis.nodes)\n",
    "\n",
    "    xplot = linspace(-1, 1, 500)\n",
    "    uplot = interpolate(xplot, u, basis)\n",
    "\n",
    "    fig1 = plot(xplot, ufunc.(xplot), label=\"u\", xguide=L\"x\", yguide=L\"u\")\n",
    "    plot!(fig1, xplot, uplot, label=L\"\\mathrm{I}(u)\")\n",
    "\n",
    "    fig2 = plot(xplot, uprim.(xplot), label=\"u'\", xguide=L\"x\", yguide=L\"u'\")\n",
    "    plot!(fig2, xplot, interpolate(xplot, basis.D*u, basis), label=L\"\\mathrm{I}(u)'\")\n",
    "\n",
    "    display(basis)\n",
    "    display(plot(fig1, fig2))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integration\n",
    "\n",
    "The nodes and weights are from [FastGaussQuadrature.jl](https://github.com/ajt60gaibb/FastGaussQuadrature.jl)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "ufunc(x) = sinpi(x)^6\n",
    "\n",
    "function compute_error(p, basis_type)\n",
    "    basis = basis_type(p)\n",
    "    u = ufunc.(basis.nodes)\n",
    "    abs(5/8 - integrate(u, basis))\n",
    "end\n",
    "\n",
    "ps = 1:23\n",
    "scatter(ps, compute_error.(ps, LobattoLegendre), label=\"Lobatto\", xguide=L\"p\", yguide=\"Error\", yaxis=:log10)\n",
    "scatter!(ps, compute_error.(ps, GaussLegendre), label=\"Gauss\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation of Orthogonal Polynomials\n",
    "\n",
    "## Legendre Polynomials\n",
    "\n",
    "Legendre poylnomials $P_p$ are evaluated as `legendre(x, p)` using the three term recursion formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "x = linspace(-1, 1, 10^3)\n",
    "fig = plot(xguide=L\"x\")\n",
    "for p in 0:5\n",
    "    plot!(fig, x, legendre.(x, p), label=\"\\$ P_$p \\$\")\n",
    "end\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gegenbauer Polynomials\n",
    "\n",
    "Gegenbauer poylnomials $C_p^{(\\alpha)}$ are evaluated as `gegenbauer(x, p, α)` using the three term recursion formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "α = 0.5\n",
    "x = linspace(-1, 1, 10^3)\n",
    "fig = plot(xguide=L\"x\")\n",
    "for p in 0:5\n",
    "    plot!(fig, x, gegenbauer.(x, p, α), label=\"\\$ C_$p^{($α)} \\$\")\n",
    "end\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jacobi Polynomials\n",
    "\n",
    "Jacobi poylnomials $P_p^{\\alpha,\\beta}$ are evaluated as `jacobi(x, p, α, β)` using the three term recursion formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "α, β = -0.5, -0.5\n",
    "#α, β = 0, 0\n",
    "x = linspace(-1, 1, 10^3)\n",
    "fig = plot(xguide=L\"x\")\n",
    "for p in 0:5\n",
    "    plot!(fig, x, jacobi.(x, p, α, β), label=\"\\$ P_$p^{$α, $β} \\$\")\n",
    "end\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hermite Polynomials\n",
    "\n",
    "Hermite poylnomials $H_p$ are evaluated as `hermite(x, p)` using the three term recursion formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "x = linspace(-2, 2, 10^3)\n",
    "fig = plot(xguide=L\"x\")\n",
    "for p in 0:4\n",
    "    plot!(fig, x, hermite.(x, p), label=\"\\$ H_$p \\$\")\n",
    "end\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hahn Polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "using LaTeXStrings, Plots; pyplot()\n",
    "\n",
    "α, β, N = -0.5, -0.5, 20\n",
    "x = linspace(0, N, 10^3)\n",
    "fig = plot(xguide=L\"x\")\n",
    "for p in 0:5\n",
    "    plot!(fig, x, hahn.(x, p, α, β, N), label=\"\\$ Q_$p(x; $α, $β, $N) \\$\")\n",
    "end\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic Computations\n",
    "\n",
    "Symbolic computations using [SymPy.jl](https://github.com/JuliaPy/SymPy.jl) and [SymEngine.jl](https://github.com/symengine/symengine) are supported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "import SymPy\n",
    "import SymEngine\n",
    "\n",
    "x_sympy = SymPy.symbols(\"x\")\n",
    "x_symengine = SymEngine.symbols(\"x\")\n",
    "\n",
    "legendre(x_sympy, 6) |> SymPy.expand |> display\n",
    "legendre(x_symengine, 6) |>  SymEngine.expand |> display\n",
    "\n",
    "GaussLegendre(2, SymPy.Sym).D |> display\n",
    "GaussLegendre(2, SymEngine.Basic).D |> display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modal Matrices\n",
    "\n",
    "The Vandermonde matrix $V$ can be computed as `legendre_vandermonde(basis)` and used to transform coeficients in a modal basis of Legendre polynomials to a nodal basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "\n",
    "p = 7\n",
    "basis = GaussLegendre(p)\n",
    "V = legendre_vandermonde(basis)\n",
    "\n",
    "# the modal derivative matrix\n",
    "Dhat = legendre_D(p)\n",
    "# they should be equal\n",
    "norm( basis.D - V * Dhat / V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the Vandermonde matrix with respect to the Jacobi polynomials is given as `jacobi_vandermonde(basis, α, β)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Revise\n",
    "using PolynomialBases\n",
    "\n",
    "p = 3\n",
    "α, β = -0.5, -0.5\n",
    "basis = GaussJacobi(p, α, β)\n",
    "V = jacobi_vandermonde(basis, α, β)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
